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Abstract. - We investigate metal- insulator transitions in the half-filled Holstein-Hubbard 
model as a function of the on-site electron-electron interaction U and the electron-phonon cou- 
pling g. We use several different numerical methods to calculate the phase diagram, the results 
of which are in excellent agreement. When the electron-electron interaction U is dominant the 
transition is to a Mott-insulator; when the electron-phonon interaction dominates, the transi- 
tion is to a localised bipolaronic state. In the former case, the transition is always found to 
be second order. This is in contrast to the transition to the bipolaronic state, which is clearly 
first order for larger values of U. We also present results for the quasiparticle weight and the 
double-occupancy as function of U and g. 



Strong correlation effects and localization can occur in metallic systems due both to strong 
electron-electron interactions and strong electron-phonon coupling and also their interplay 
There are many systems with strongly correlated electrons where there is also a strong coupling 
to the lattice and lattice modes, for example V2O3 [1,2], manganites [3] or fullerides [4]. 
The strong electron-electron interactions can be described by the Hubbard model, where the 
transition to a Mott-Hubbard insulator has been extensively investigated [5-9] . The Holstein 
model has been used to examine localization to a polaronic or bipolaronic insulator due to 
electron-phonon interactions [10-13] . The interplay between these two types of localization can 
be investigated using the Holstein-Hubbard model which includes both types of interaction: 

ka » i i <T 

where U describes the electron-electron interaction within a semielliptic band of dispersion 
e(k) and band width W = 4, g the electron-phonon coupling and ujq = 0.2 is the frequency 
of a local Einstein-phonon mode. In the large-wo limit (cjq ^> W), the model can be mapped 
onto an effective Hubbard model with U c s = U — 2g 2 /ujo [14, 15]. 

Here we investigate the phase diagram in the particle-hole symmetric case. We neglect 
phases of long-range order and concentrate on metal to insulator transitions driven by both 
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Fig. 1 - The calculated phase diagram of the Hubbard-Holstein model at T = 0. The solid lines 
represent the DIA2 results. They are in excellent agreement with the DIA4, NRG, and ED results, 
as indicated by the crosses, open circles, and filled circles respectively. The dashed line represents the 
locus of the points with U c g = U — 2g 2 /u>o = which becomes the phase boundary for U > 6.4. 



electron-electron and elcctron-phonon interactions [16]. In general, no exact solution of this 
model is known. We use a number of approximation schemes which lead to a consistent picture 
of the physics of these transitions. Two of these approaches are based on the dynamical mean- 
field theory [7], where the lattice model is mapped onto an effective impurity model with a 
self-consistency constraint. We apply two different methods to solve the effective impurity 
model, the exact diagonalization method (ED) [17] and the numerical renormalization group 
(NRG) [18, 19]. In the ED method a restricted basis set is used to describe the bath of the 
impurity model, which can then be solved exactly. The ED calculations presented here were 
performed on an 8-site cluster. In the NRG approach, the impurity model is solved by an 
iterative diagonalization scheme which has the advantage that it can probe arbitrarily low 
energy scales. In the NRG calculations up to 1200 states were retained with A = 1.8 and up 
to 36 phonon states. We also use the dynamical impurity approximation (DIA), introduced 
recently by Potthoff [20,21]. This method has the advantage of being a variational approach 
for the grand potential. Its main limitation is the restricted set of trial selfenergies that one 
can handle in practice. We use trial selfenergies corresponding to impurity models with one 
(two-site DIA - DIA2) and three bath sites (four-site DIA - DIA4). 

In figure^ we plot the g vs. U phase diagram as obtained by the methods described above. 
There is excellent agreement between the results of all four methods. The critical coupling for 
the Mott transition U c m ~ 5.85 on the g = axis corresponds to the results already known 
for the Hubbard model [8,9], and g c w 0.39 on the U = axis for the Holstein model [12,13]. 
Within the DMFT for U < U c m there is a range of values of U where the insulating solution 
coexists with the metallic one (" hysteresis" ) , but for T — the metallic solution is always the 
physically relevant one. For the Holstein model a significant coexistence region exists only 
for large loq [15]. The dashed line (polaronic line) is the locus of the points with U e g = 0; 
above (below) this line, U c s < (> 0). The metal-insulator transition in the region U c h > 
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Fig. 2 - The quasiparticle weight z as a function of the electron-electron repulsion U for two values 
of the electron-phonon coupling: g — 0.3 in the left figure and g = 0.45 in the right. 

is very similar to the Mott-Hubbard transition as found in the pure Hubbard model, where 
double-occupancy d = {n^n^) is almost completely suppressed. 

The phase boundary to the Mott-Hubbard insulator is largely independent of the electron- 
phonon coupling until the polaronic line is reached. This is due to the fact that in this region 
the charge fluctuations which couple to the phonons are strongly suppressed. On the other 
hand, the transition in the other region (U c s < 0) depends strongly on both U and g. The 
insulating phase here is a bipolaronic one with enhanced on-site double-occupancy. 

To probe the physics of these transitions in more detail, we consider the variation of the 
quasiparticle weight z both in scans with a fixed g or a fixed U. In figurc|5|we look at the cases 
g = 0.3 and g — 0.45. The results of all four methods show the same general trend, although 
the two-site DIA appears to systematically overestimate the values of z in general but gives 
almost correct values for the critical couplings. For g = 0.3 we start in a metallic state for 
[7 = with a value of z which is already rather less than one, due to the electron-phonon 
interaction [12]. Upon increasing U, z initially increases until U e s rj and then steadily 
decreases to z = at the critical value of U = U c m- For g = 0.45 and U = we start in 
the insulating phase. At a critical value U = U c b ~ 0.9 there is a sudden onset of a metallic 
solution which suggests that this transition might be first order. As in the case for g = 0.3, the 
quasiparticle weight increases until U e s rs and decreases to the transition at U = U c m ~ 6.0. 
In both cases, z vanishes continuously as U approaches U c m from below. 

The complementary scans of z with fixed U and variable g are shown in figure 01 There 
is little change of z in the region with U c g > as the repulsive electron-electron interaction 
U effectively suppresses charge fluctuations which would otherwise couple to the phonons. 
This suppression can also be seen in the results for the double-occupancy d — (njni), which 
are plotted in figure |3 Once U e g is less than zero, there is a rather rapid decrease in z 
and an apparent jump at the critical value of g indicating a first order transition. There is 
a corresponding sudden increase of the double-occupancy d. For larger values of U, z even 
increases until a point very close to the transition, and the jump becomes more pronounced 
as can be seen for U = 5.0 (lower left figure EJ. Above the transition the value of d s=s 0.5. 

With numerical methods it is difficult to exclude the possibility that the transition is 
continuous but very sharp. However, with the DIA one calculates the grand potential £1 
which enables one to address the order of the transition directly. In figure 01 we plot the DIA2 
results for as a function of the variation parameter V for g = 0.60 and values of U close 
to the transitions at U c m (left plot) and U c b (right plot). For the transition at U c m we see 
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Fig. 3 - The quasiparticle weight z as a function of the electron-phonon coupling g for three values of 
the electron-electron interaction: U = 1.0 in the upper left, U = 3.0 in the upper right and U = 5.0 
in the lower left figure. The lower right figure shows the double-occupancy d = (n^n\) as a function 
of g for all three values of U. The lines correspond to the DIA2 results and the open (filled) circles 
to the NRG (ED) results. 



that there is a global minimum which shifts continuously to V — as U approaches U c m- In 
the DIA2, V — corresponds to the insulating solution. These results are very similar to the 
results obtained by Potthoff for the pure Hubbard model [21]. 

In the right-hand figure the same data is plotted for the transition at U c b ■ Starting with the 
largest value of U = 3.26 the global minimum is at finite V w 0.56 corresponding to a metallic 
state. Decreasing U, the minimum becomes shallower but shifts only slightly to smaller values 
V w 0.54. The minimum at V = becomes the global minimum for U < U c b = 3.225. This 
is clear evidence for a first order transition, at least within the DIA2. Coexistence of metallic 
and insulating solutions can be found there, as is corroborated by the other methods. 

In figure |5] we plot the grand potential as a function of g as calculated with the DIA2 and 
DIA4 approaches for three values of U = 1.0,3.0 and 5.0. For the largest value of U = 5.0, 
the discontinuity of the gradient of O with respect to U is clearly seen. As one reduces U, 
the kink becomes progressively smaller. For DIA2 calculations, similar to those presented in 
figure a kink persists down to U = and the phase transition is always first order. For 
the corresponding DIA4 calculations, however, the situation is not quite so clear. In this 
small-[/ range, a jump in z is also difficult to identify (see figure [3J). The evidence from the 
calculations apart from the DIA2, indicates that the transition is probably for all practical 
purposes continuous for U < 3. In contrast, the transitions along the line U c m are always 
continuous as a function of U. 
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Fig. 4 - The grand potential Q as a function of the variational parameter V of the DIA2 method at 
g = 0.60. The left figure displays a second order transition from a metal to a Mott-Hubbard insulator. 
The minimum shifts smoothly to zero as U increases from U = 5.6, the lowest curve, in steps of 0.2 
to U = 7.0 for the uppermost curve. The right figure shows a first order transition from a metal to 
a bipolaronic state. There is a jump of the position of the minimum as U decreases from U = 3.26, 
the lowest curve, in steps of —0.01 to U = 3.20 for the uppermost curve. 



All methods used in our investigation lead to a consistent picture of the metal-insulator 
transitions in the Holstein-Hubbard model. The phase diagram is composed of three regions, 
a metallic region and two qualitatively different insulating regions. The boundary between the 
metallic and the Mott-insulating region is described by a line of continuous transitions, and is 
only weakly dependent on the electron-phonon coupling. Our evidence is that the transition 
from the metal to a bipolaron insulator state, which is characterised by an enhanced on-site 
double-occupancy, is of first order for larger values of U. However, the discontinuity becomes 




Fig. 5 - The grand potential £1 as function of g for three different values of U , as calculated with the 
DIA2 (big open symbols) and DIA4 (small dots) methods. Here, Q,2{V = 0) denotes the value of the 
grand potential for V = in the DIA2 method. 
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progressively smaller as U is decreased and has virtually disappeared for the range U < 3. 

The behaviour near the boundary between the two insulating phases is difficult to access 
by all of our methods except within the DIA2, where the results indicate that there is a 
first order transition between them. Detailed results and discussions of dynamical response 
functions at T — for this model have been calculated and will be prepared for publication 
in due course. 

* * * 

We wish to thank the EPSRC (Grant GR/S18571/01) for financial support. One of us (YO) 
was supported by the Grant-in-Aid for Scientific Research from the Ministry of Education, 
Culture, Sports, Science and Technology. We also thank M. Aichhorn, R. Bulla and M. 
Potthoff for helpful discussions. 

REFERENCES 

[1] P. Majumdar and H. R. Krishnamurthy, Phys. Rev. Lett. 73(11), 1525 1994. 

[2] M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. 70, 1039 1998. 

[3] A. J. Millis, P. B. Littlewood, and B. I. Shraiman, Phys. Rev. Lett. 74(25), 5144 1995. 

[4] O. Gunnarson, Rev. Mod. Phys. 69(2), 575 1997. 

[5] J. Hubbard, Proc. R. Soc. London, Ser. A 281, 401 1964. 

[6] M. Jarrell, Phys. Rev. Lett. 69(1), 168 1992. 

[7] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68(1), 13 1996. 

[8] G. Moeller, Q. Si, G. Kotliar, M. Rozenberg, and D. S. Fisher, Phys. Rev. Lett. 74, 2082 1995. 

[9] R. Bulla, Phys. Rev. Lett. 83, 136 1999. 

[10] T. Holstein, Ann. Phys. 8, 325 1959. 

[11] P. Benedetti and R. Zeyher, Phys. Rev. B 58(21), 14320 1998. 

[12] D. Meyer, A. C. Hewson, and R. Bulla, Phys. Rev. Lett. 89(19), 196401 2002. 

[13] M. Capone and S. Ciuchi, Phys. Rev. Lett. 91(18), 186405 2003. 

[14] J. K. Freericks, Phys. Rev. B 48(6), 3881 1993. 

[15] D. Meyer and A. C. Hewson, Acta Physica Polonica B 34, 769 2003. 

[16] Recently this problem has also been considered using the NRG approach by G. S. Jeon, T.-H. 

Park, J. H. Han, H. C. Lee, and H.-Y. Choi, cond-mat/0312390. 

[17] M. Caffarel and W. Krauth, Phys. Rev. Lett. 72(10), 1545 1994. 

[18] H. R. Krishna-murthy, J. W. Wilkins, and K. G. Wilson, Phys. Rev. B 21(3), 1003 1980. 

[19] R. Bulla, A. C. Hewson, and T. Pruschke, J. Phys.: Condens. Matter 10, 8365 1998. 

[20] M. Potthoff, Eur. Phys. J. B 32, 429 2003. 

[21] M. Potthoff, Eur. Phys. J. B 36, 335 2003. 



